//Port of JAMPACK's householder implementation. 
using System;
using System.Collections.Generic;
using System.Text;

using dnAnalytics.Exceptions;
using dnAnalytics.LinearAlgebra;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal static class Householder
    {
        public static double Norm(double[] vector)
        {
            double ret = 0;
            for (int i = 0; i < vector.Length; ++i)
            {
                ret += vector[i] * vector[i];
            }
            return System.Math.Sqrt(ret);
        }

        public static double[] GenerateColumn(double[] A, int m, int r1, int r2, int c)
        {
            int ru = r2 - r1 + 1;
            int tmp = c * m;
            int index = tmp + r1;
            double[] u = new double[ru];

            for (int i = r1; i <= r2; i++)
            {
                u[i - r1] = A[tmp + i];
                A[tmp + i] = 0.0;
            }


            double norm = Norm(u);

            if (r1 == r2 || norm == 0)
            {
                A[index] = -u[0];// A.ValueAt(r1, c, -u[0]);
                u[0] = System.Math.Sqrt(2.0);
                return u;
            }

            double scale = 1.0 / norm;
            if (u[0] < 0.0)
            {
                scale *= -1.0;
            }

            A[index] = -1.0 / scale;// A.ValueAt(r1, c, (-1.0 / scale));

            for (int i = 0; i < ru; i++)
            {
                u[i] *= scale;
            }

            u[0] += 1.0;
            double s = System.Math.Sqrt(1.0 / u[0]);

            for (int i = 0; i < ru; i++)
            {
                u[i] *= s;
            }
            return u;
        }

        public static void UA(double[] u, double[] A, int m, int r1, int r2, int c1, int c2, double[] v)
        {
            if (r2 < r1 || c2 < c1)
            {
                return;
            }

            if (r2 - r1 + 1 > u.Length)
            {
                throw new NotConformableException();
            }

            if (c2 - c1 + 1 > v.Length)
            {
                throw new NotConformableException();
            }

            for (int j = c1; j <= c2; j++)
            {
                v[j - c1] = 0.0;
            }

            for (int i = r1; i <= r2; i++)
            {
                for (int j = c1; j <= c2; j++)
                {
                    v[j - c1] = v[j - c1] + u[i - r1] * A[j * m + i];
                }
            }

            for (int i = r1; i <= r2; i++)
            {
                for (int j = c1; j <= c2; j++)
                {
                    A[j * m + i] = A[j * m + i] - u[i - r1] * v[j - c1];
                }
            }
        }

        public static void UA(double[] u, double[] A, int m, int r1, int r2, int c1, int c2)
        {
            if (c1 > c2)
            {
                return;
            }
            UA(u, A, m, r1, r2, c1, c2, new double[c2 - c1 + 1]);
        }

    }
}
